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Implicit  high  order  strong  stability  preserving  Runge-  Kutta  time 

DISCRETIZATIONS 

AFOSR  GRANT  NUMBER  FA9550-06-1-0255  Sigal  Gottlieb 
Mathematics  Department 
University  of  Massachusetts  -  Dartmouth 


Abstract 

We  investigated  diagonally  split  Rmige-Kutta  (DSRK)  methods  to  identify  and  test 
unconditionally  strong  stability  preserving  (SSP)  methods,  and  implicit  SSP  time- 
stopping  methods  to  find  methods  with  a  large  SSP  eoeffieient.  We  found  that  DSRK 
methods  which  are  unconditionally  SSP  reduce  to  first  order  for  the  stepsizes  of 
interest,  and  introduced  an  analysis  which  explains  this  phenomenon  and  shows  that  it 
is  unavoidable.  We  found  optima;  implicit  SSP  Runge  Kutta  methods  up  to  order  six 
(which  is  the  maximal  possible  order  for  these  methods)  and  eleven  stages,  and  found 
that  the  eifcctive  SSP  coefficient  can  be  no  more  than  two,  making  these  methods 
not  competitive  with  explicit  methods  for  most  applications,  but  useful  in  a  carefully 
chosen  subset  of  problems.  We  now  have  a  complete  analysis  of  implicit  SSP  Runge 
Kutta  methods  and  demonstrations  of  the  need  for  the  SSP  property  in  solutions  of 
hyperbolic  PDFs  with  shocks. 


1  Summary  of  Aims  and  Results 


Strong  stability  preserving  (SSP)  high  order  time  discretizations  were  developed  to 
ensure  nonlinear  stability  properties  necessary  in  the  numerical  solution  of  hyperbolic 
partial  differential  equations  with  discontinuous  solutions.  SSP  methods  preserve  the 
strong  stability  properties  -  in  any  norm,  serninorni  or  convex  functional  -  of  the  spa¬ 
tial  discretization  coupled  with  first  order  Euler  time  stepping,  when  the  timestep  is 
suitably  restricted.  Explicit  strong  stability  preserving  (SSP)  Runge-Kutta  methods 
([17],  [18], [19],  [20],  [4], [5],  [6])  have  been  successfully  used  with  a  wide  range  of  spatial 
discretizations,  including  spectral,  discontinuous  Galerkiii,  and  weighted  essentially 
non-oscillatory  (WENO)  methods.  These  high  order  methods  preserve  any  nonlinear 
stability  properties  satisfied  by  the  spatial  discretization  coupled  with  the  forward 
Euler  time-stepping.  However,  all  general  linear  methods  suffer  from  a  SSP  time-step 
restriction.  This  motivates  the  search  for  high  order  implicit  time-stepping  methods 
with  SSP  properties  and  a  large  allowable  time-step,  which  is  the  overarching  goal  of 
this  project. 

The  connections  between  the  SSP  property  and  the  theory  of  contractivity  have 
provided  efficient  tools  for  the  study  of  SSP  multistep  and  Runge-Kutta  methods, 
which  we  utilized  in  the  search  for  optimal  implicit  SSP  Runge-Kutta  methods.  Fur¬ 
thermore,  contractivity  theory  allowed  us  to  determine  order  barriers  on  SSP  methods. 
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to  establish  l:)oiiiids  on  the  SSP  coefficient,  and  to  conclude  that  the  SSP  coefficient 
is  not  only  sufficient  but  necessary  for  strong  stability  preservation  in  an  arbitrary 
norm  for  an  arbitrary  semi-discretization  that  satisfies  a  strong  stability  condition 
under  forward  Euler  integration. 

This  work  described  below  was  performed  in  collaboration  with  David  Ketcheson, 
a  doctoral  candidate  at  University  of  Washington  in  Seattle  (advised  by  Dr,  R.  LeV- 
eque),  and  Colin  Macdonald,  then  a  doctoral  candidate  at  Simon  Fraser  University 
(advised  by  Dr.  S.  Ruuth)  and  now  a  postdoctoral  fellow  at  UCLA  (working  with 
Dr.  S.  Osher). 

The  aims  of  AFOSR  grant  number  FA9550-06- 1-0255  were  to 


•  Use  results  and  formulations  from  contractivity  and  monotoiiicity  theory  to  find 
optimal  class  of  higher  order  implicit  SSP  Ruiige-Kutta  methods. 

•  Find  higher  order  implicit  diagonally  split  Runge-Kutta  methods  (DSRK)  which 
are  SSP  methods  with  no  stepsize  restriction. 

•  Test  the  optimal  implicit  SSP  Runge-Kutta  methods  for  use  with  flux-implicit 
WENO  spatial  discretizations. 

•  Test  the  DSRK  with  no  time-step  restriction  with  spectral  and  WENO  spatial 
discretizations. 


In  the  grant  period  we  have  gone  further  than  we  proposed  or  anticipated.  The 
following  are  the  accoinplisliments  under  this  grant: 


1.  We  conducted  a  thorough  numerical  study  of  second  and  third  order  diagonally 
split  Runge-Kutta  methods  on  a  variety  of  problems.  These  methods  have 
proved  disappointing,  due  to  severe  reduction  of  order  which  renders  them  no 
better  than  backward  Euler,  which  is  unconditionally  SSP  [16].  We  analyzed 
the  cause  of  this  order  reduction  and  found  a  way  to  avoid  it,  however  this 
rendcTs  the  SSP  coefficient  as  small  as  for  implicit  Runge-Kutta  methods. 

2.  The  connections  between  the  SSP  property  and  the  theory  of  contractivity 
optimal  have  provided  efficient  tools  for  the  study  of  SSP  multistep  and  Runge- 
Kutta  methods.  Methods  of  these  types  have  been  thoroughly  investigated,  and 
their  development  seems  to  be  essentially  complete.  Furthermore,  contractiv¬ 
ity  theory  allowed  us  to  determine  order  barriers  on  SSP  methods,  to  establish 
bounds  on  the  SSP  coefficient,  and  to  conclude  that  the  SSP  coefficient  is  not 
only  sufficient  but  necessary  for  strong  stability  preservation  in  an  arbitrary 
norm  for  an  arbitrary  semi-discretization  that  satisfies  a  strong  stability  condi¬ 
tion  under  forward  Euler  integration. 
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3.  We  found  optimal  implicit  SSP  Runge-Kutta  methods  up  to  order  six  and  up 
to  eleven  stages.  These  methods  are  diagonally  implicit  or  singly  diagonally 
implicit  and  have  sparse,  efficient  representations.  Furthermore,  the  implicit 
solutions  at  each  stage  of  a  SSP  Runge-Kutta  method  have  provable  existence 
and  uniqueness  properties. 

4.  Our  work  demonstrated  that  implicit  SSP  methods  are  unlikely  to  bo  efficic'iit 

enough  to  out-perform  the  explicit  methods.  We  define  the  effective  SSP  co- 
efficient  of  a  method  ^  to  normalize  the  step-size  coefficient  c  relative 

to  the  number  of  stages  m  in  a  method.  The  very  restrictive  bound  Ce//  <  2 
has  been  proven  for  implicit  multistep  methods  [15,  10]  and  conjectured  for  im¬ 
plicit  Runge-Kutta  methods  [12].  In  contrast,  explicit  methods  have  a  bound 
^ef f  —  I* 

5.  In  the  wider  class  of  explicit  general  linear  methods  (which  includes  both  Ruiige 
Kutta  and  multistep  methods  as  a  subset)  the  bound  Ce//  <  I  was  proved  [7]. 

6.  Although  the  focus  of  this  grant  was  implicit  Runge-Kutta  methods,  the  tools 
developed  for  this  grant  allowed  David  Ketcheson  to  independently  perform  a 
more  thorough  study  of  explicit  low-storage  Runge-Kutta  methods  [13]  as  well 
as  implicit  and  explicit  multi-step  methods  [7].  We  found  that  the  SSP  Runge- 
Kutta  methods  tend  to  have  a  variety  of  nice  properties,  such  as  small  error 
constants  and  large  regions  of  absolute  stability. 

7.  We  showed  that  spectral  deferred  correction  methods  can  be  written  as  Runge- 
Kutta  method,  and  are  thus  amenable  to  the  techniques  for  efficient  optimiza¬ 
tion  found  using  the  connections  to  contractivity  theory.  Using  these  connec¬ 
tions,  we  also  conclude  that  these  methods  suffer  from  the  same  order  barriers 
and  bounds  on  the  SSP  coefficient. 

8.  David  Ketcheson  further  studied  the  SSP  properties  of  the  Runge-KuttaChebyshev 
methods.  Verwers  second  order  methods  all  have  negative  Butcher  coe?cients, 

so  they  are  not  SSP  under  any  positive  timestep.  We  have  found  first  and  second 
order  SSP  methods  up  to  10  stages  that  have  the  theoretically  optimal  time- 
step.  These  are  promising  for  fully  explicit  integration  of  convection-di?usioii 
equations  without  operator  splitting.  Unlike  IMEX,  exponential  di?ereiicing, 
etc.,  they  apply  the  same  integration  method  to  the  sti?  and  non-sti?  parts) 


1.1  Publications: 

Publications  resulting  from  this  grant  are; 

1.  “A  numerical  study  of  diagonally  split  Runge-Kutta  methods  for  PDEs  with 
discontinuities’’  by  C.B.  Macdonald,  S.  Gottlieb,  and  S.  Ruuth.  Journal  of 
Scientific  Computing,  36(1):89-112,  (2008). 
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2.  “Optimal  implicit  strong  stability  preserving  Runge-Kutta  methods’'  by  D. 
Ketclieson,  C.  Macdonald,  and  S.  Gottlieb.  Applied  Numerical  Mathematics 
(to  appear). 

3.  “Highly  E?cieiit  Strong  Stability  Preserving  Runge-Kutta  Alethods  with  Low 
Storage  Implementations”  by  D.  Ketclieson.  SIAM  Journal  on  Scieiitic  Com¬ 
puting,  30  (4):  2113-2136  (2008).  Winner  of  the  SIAM  student  paper  prize. 

4.  “Computation  of  optimal  monotonicity  preserving  general  linear  methods”  by 
David  L  Ketclieson.  Math,  of  Comp.  (2008) 

5.  “High  Order  Strong  Stability  Preserving  Time  Discretizations”  by  S.  Gottlieb. 
D.I.  Ketcheson  and  C.-W.  Shu.  Journal  of  Scientific  Computing  38:251-289 
(2009). 


1.2  Dissemination 


Other  dissemination  efforts  related  to  this  grant: 


1.  We  set  up  a  web-site  devoted  to  SSP  methods,  to  collect  all  the  latest  results  and 
most  useful  information  about  strong  stability  preserving  time  discretizations, 
http:/ /www.cfm. brown.edu/people/sg/ssp.htinl 

2.  We  organized  a  minisymposiuin  at  the  2006  annual  SIAM  conference  which 
brought  together  Rong  Wang  (who  presented  his  joint  work  with  Ray  Spiteri), 
luma  Higueras,  Steven  Ruuth  and  his  student  Colin  Macdonald.  This  niinisyin- 
posia  led  to  productive  discussions  with  Adrian  Sandu  and  his  student  on  the 
topic  of  SSP  multirate  time-stepping. 

(a)  Positivity  and  Monotonicity  for  IMEX  Methods  by  Inmaculada  Higueras, 
Uiiiversidad  Pblica  de  Navarra,  Spain. 

(b)  Variable  Step-Size  IMplicit-EXplicit  Linear  Multistep  Methods  by  Steve 
Ruuth,  Simon  Fraser  University,  Canada;  Dong  Wang,  University  of  Illi¬ 
nois  at  Urbana-Champaign. 

(c)  In  Search  of  Implicit  High-Order  Strong  Stability  Preserving  Methods  with 
Relaxed  Time-Step  Restrictions  Sigsl  Gottlieb,  University  of  Massachusetts: 
Colin  Macdonald,  Simon  Fraser  University;  Steve  Ruuth,  Simon  Fraser 
University,  Canada. 

(d)  Comments  on  Linear  Instability  of  Time  Integration  Methods  with  the 
Fifth- Order  WENO  Spatial  Discretization  Raymond  J.  Spiteri  and  Rong 
Wang,  University  of  Saskatchewan,  Canada. 

3.  We  have  organized  a  minisymposium  which  will  take  place  at  the  2008  SIAM 
annual  meeting,  which  will  feature  the  following: 
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(a)  Strong  Stability  Preserving  Time- Stepping  Methods  by  Sigal  Gottlieb. 
University  of  Massachusetts,  Dartmouth;  David  Ketcheson,  University  of 
Washington;  Colin  Macdonald,  Simon  Fraser  University 

(b)  Optimal  Explicit  and  Implicit  SSP  Runge-Kutta  Methods  by  David  I. 
Ketcheson,  University  of  Washington;  Colin  Macdonald,  Simon  Fraser 
University;  Sigal  Gottlieb,  University  of  Massachusetts,  Dartmouth; 

(c)  Practical  considerations  for  IMEX  SSP  Runge-Kutta  methods  by  Inmac- 
ulada  Higueras,  Universidad  Pblica  de  Navarra,  Spain;  Teo  Roldan. 

(d)  Generalizations  of  Positivity  and  Strong  Stability  Preservation  by  Zoltan 
Horvath,  Szcheiiyi  Istvan  University,  Gyr,  Hungary. 

(e)  High  Order  Discretizations  of  Kinetic  Equations  by  Lorenzo  Pareschi. 
University  of  Ferrara,  Italy. 

(f)  Multirate  SSP  Methods  for  Hyperbolic  PDEs  by  Emil  Constantinescu 
and  Adrian  Sandu,  Virginia  Polytechnic  Institute  &:  State  University. 

(g)  Do  We  Know  WENO?  by  Raymond  J.  Spiteri  and  Rong  Wang,  Uni¬ 
versity  of  Saskatchewan,  Canada. 

(h)  Stage-exceeding  Order  SSP  Time-stepping  for  Runge-Kutta  Discontinuous 
Galerkin  Methods  by  Clint  Dawson,  University  of  Texas,  Austin;  Ethan 
Kubatko,  University  of  Texas  at  Austin. 

4.  Seminar  presentation  "Strong  Stability  Preserving  time  discretizations  with  op¬ 
timal  time-step  restrictions”  at  UMass  Amherst  on  October  30,  2007. 

5.  Workshop  presentation  "Strong  Stability  Preserving  Time  Discretizations”  at 
the  Statistical  and  Applied  Mathematical  Sciences  Institute’s  (SAMSI)  2007- 
2008  Program  on  Random  Media  Interface  Problems  Workshop  in  North  Car¬ 
olina  on  November  15,  2007. 

6.  Seminar  presentation  "On  Strong  Stability  Preserving  Runge-Kutta  and  Multi- 
step  Time  Discretizations"  at  the  (NYU)  Courant  Institute’s  Numerical  Anal¬ 
ysis  and  Scientific  Computing  seminar  on  November  30,  2007. 

7.  Seminar  presentation  “Time  stepping  methods  for  numerical  solution  of  hy¬ 
perbolic  PDEs  with  shocks”  in  MIT’s  Mathematics  Department’s  Numerical 
Methods  for  Partial  Differential  Equations  seminar  on  November  12,  2008. 

8.  Seminar  presentation  “Time  stepping  methods  for  numerical  solution  of  hyper¬ 
bolic  PDEs  with  shocks"  Mathematics  Department  Colloquium  in  The  Univer¬ 
sity  of  Connecticut  -  Storrs  on  November  13,  2008. 

9.  Book  contract  with  for  World  Scientific  Publishing  for  a  monograph  on  SSP 
time-discretization  methods  (together  with  C.-W.  Shu). 
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Tliese  iiiiiiisyiiiposia,  seminars,  workshops,  and  website  have  caused  the  topic  of  time- 
stepping  to  be  more  widely  discussed  and  studied,  and  has  inspired  collaborations  and 
other  research  on  the  topic. 


2  Detailed  Progress  By  Year: 


Year  1:  We  studied  the  class  of  diagonally  split  Runge-Kutta  methods  to  find  high 
order,  unconditionally  SSP  methods.  Diagonally  split  Runge-Kutta  (DSRK)  ([1, 
2,  8,  9])  time  discretization  methods  are  a  class  of  implicit  time-stepping  schemes 
which  offer  both  (formal)  high-order  convergence  and  a  form  of  nonlinear  stability 
known  as  unconditional  contractivity.  This  combination  is  not  possible  within  tlu' 
classes  of  Runge-  Kutta  or  linear  multistep  methods  and  therefore  appears  promising 
for  the  strong  stability  preserving  (SSP)  time-stepping  community  which  is  generally 
concerned  with  computing  oscillation-free  numerical  solutions  of  PDEs. 

We  conducted  a  thorough  numerical  study  of  second  and  third  order  diagonally 
split  Rurige  Kutta  methods  on  a  variety  of  of  archetypal  test  cases  including  linear 
advection.  Burgers’  equation,  a  diffusion  equation  with  discontinuous  initial  data, 
and  the  Black-Scholes  equation.  The  numerical  tests  verified  the  asymptotic  order 
of  the  schemes  as  well  as  the  unconditional  contractivity  property.  However,  in  ev¬ 
ery  numerical  experiment,  diagonally  split  Runge-Kutta  methods  suffer  from  order 
reduction  at  large  step-sizes.  Indeed,  for  time-steps  larger  than  those  typically  cho¬ 
sen  for  explicit  methods,  these  diagonally  split  Runge-Kutta  methods  behave  like 
first-order  implicit  methods.  In  every  numerical  experiment,  the  unconditionally  con¬ 
tractive  diagonally  split  Runge-Kutta  methods  were  out-performed  by  the  first-order 
backward  Euler  scheme  when  At  >  2AtFEi  by  explicit  Runge-Kutta  methods  or 
Crank-Nieolson  when  At  <  2AtfE-  At  larger  time-steps,  the  unconditionally  con¬ 
tractive  diagonally  split  Runge-Kutta  schemes  are  strong  stability  preserving  (SSP) 
but  suffer  from  order  reduction,  making  backward  Euler  a  better  choice.  At  small 
step-sizes,  Crank-Nicolson  and  explicit  SSP  Runge-Kutta  methods  are  SSP,  and  pro¬ 
duce  far  more  accurate  results  at  a  smaller  computational  cost.  Indeed,  for  time-steps 
larger  than  those  typically  chosen  for  explicit  methods,  these  DSRK  methods  behave 
like  first-order  implicit  methods.  This  is  unfortunate,  because  it  is  precisely  to  al¬ 
low  a  large  time-step  that  we  choose  to  use  implicit  methods.  We  studied  this  order 
reduction  phenomenon  analytically,  and  showed  that  higher  stage  order  of  the  un¬ 
derlying  Runge-Kutta  schemes  was  insufficient  to  avoid  order  reduction.  We  then 
derived  DSRK  stage  order  conditions  and  constructed  DSRK  schemes  with  higher 
stage  which  do  not  suffer  from  order  reduction.  However,  because  of  the  high  stage 
order,  these  schemes  cannot  be  unconditionally  contractive,  and  the  resulting  SSP 
coefficient  are  comparable  to  implicit  Runge-Kutta  [16]. 

Year  2:  In  the  second  year  of  the  project  we  surveyed  the  literature  on  contrac- 
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tive  methods  and  extracted  results  which  are  applicable  to  SSP  methods,  identified 
efficient  techniques  to  find  the  radius  of  absolute  monotonicity,  and  found  optimal 
implicit  SSP  Runge-Kutta  methods  of  order  up  to  six. 

Using  the  results  from  coiitractivity  theory,  we  were  able  to  identify  the  following 
order  barriers  and  bounds  on  the  SSP  coefficient  of  Runge-Kutta,  multistep,  and 
general  linear  methods: 


•  Runge-Kutta  Methods 

1.  An  SSP  Runge-Kutta  method  with  can  be  no  more  than  fourth  order 
accurate  if  it  is  explicit  and  no  more  than  sixth  order  accurate  if  it  is 
implicit  [14]. 

2.  Implicit  Runge-Kutta  methods  that  are  unconditionally  SSP  must  have 
order  at  most  one.  This  result  is  in  contreist  with  linear  stability  and 
stability,  where  some  high^order  implicit  methods  (i.e.,  the  A-stable  meth¬ 
ods  and  the  algebraically  stable  methods,  respectively)  are  unconditionally 
stable. 

3.  The  implicit  SSP  Runge-Kutta  of  order  p  >  1  have  an  SSP  coefficient  that 
is  not  dramatically  larger  than  those  for  explicit  methods  [15,  3,  12]. 

4.  Any  SSP  method  must  have  stage  order  p  <  2,  and  explicit  Runge  Kutta 
method  must  have  stage  order  p  <\.  The  stage  order  p  is  a  lower  bound  on 
the  order  of  convergence  when  a  method  is  applied  to  arbitrarily  stiff  prob¬ 
lems.  Low  stage  order  may  lead  to  oi'der  reduction,  i.e.  slow  convergence, 
when  computing  solutions  of  stiff  ODEs. 

5.  All  m-stage  diagonally  implicit  methods  have  order  at  most  rri  +  1. 

6.  All  SSP  m-stage  singly  diagonally  implicit  methods  have  order  at  most 
m  +  1 . 

7.  SSP  singly  diagonally  implicit  methods,  which  are  both  singly  implicit 
and  diagonally  implicit,  have  the  same  order  barrier  (p  <  4)  as  explicit 
methods. 

•  Multistep  Methods 

1.  For  explicit  s-step  methods  of  order  p,  the  SSP  coefficient  is  bounded  by 
c  <  ^  for  s  >  1 

2.  for  implicit  methods  of  order  p  >  1,  the  SSP  coefficient  is  bounded  by 

c  <  2. 

3.  While  there  appears  to  be  no  limit  to  the  order  of  accuracy  of  SSP  linear 
multistep  methods,  high  order  accurate  methods  of  this  type  are  subject 
to  very  small  timestep  restrictions  and  require  very  many  steps. 
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•  General  Linear  Methods 

Any  explicit  m-stage,  .s-step  general  linear  method  of  order  p,  has  SSP  coefficient 
bounded  by  the  number  of  its  stages,  c  <  m. 


Although  no  unconditionally  SSP  method  can  have  order  greater  than  one  [21],  we 
explored  the  possibility  that  implicit  methods  may  have  SSP  coefficients  signifi(*antly 
larger  than  those  of  explicit  methods  with  the  same  order  and  number  of  stages. 
The  question  we  wished  to  answer  was  whether  the  allowable  step-size  can  be  large 
enough  to  offset  the  extra  computational  effort  required  in  the  implicit  solution  of 
the  resulting  system  at  each  iteration. 

Using  the  efficient  formulation  of  the  problem  of  finding  the  radius  of  coiitractivity 
of  a  method,  it  was  possible  to  use  MATLAB  to  perform  a  search  for  optimal  implicit 
SSP  Rungc-Kutta  methods.  These  results  gave  us  optimal  methods  of  order  up  to 
six,  which  is  the  maximal  order  for  implicit  SSP  Runge-Kutta  methods.  In  fact,  only 
existence  of  methods  of  order  up  to  five  was  previously  established  [14].  Our  search 
successfully  found  methods  of  order  six,  establishing  that  this  is  indeed  possible  and 
that  the  order  barrier  is  sharp. 

Recently,  Ferracina  and  Spijker  investigated  optimal  singly  diagonally  implicit 
Runge-Kutta  methods  [3].  They  showed  that  such  methods  have  order  at  most  four, 
and  found  optimal  methods  (by  numerical  optimization)  of  up  to  order  four  and  up 
to  eight  stages.  They  also  conjectured  the  form  of  optimal  second  and  third  order 
methods  with  any  number  of  stages.  Using  numerical  optimization  techniques,  we 
performed  an  extensive  search  among  the  much  larger  class  of  fully  implicit  SSP 
Runge-Kutta  methods  [12].  Remarkably,  searching  among  the  cla^s  of  fully  implicit 
methods,  the  optimal  methods  of  second  and  third  order  were  found  to  be  singly 
diagonally  implicit;  in  fact,  they  were  the  very  methods  found  already  in  [3].  The 
optimal  methods  of  fourth  through  sixth  order  were  found  to  be  diagonally  implicit. 
Many  of  these  implicit  methods  have  representations  that  allow  for  very  efficient 
implementation  in  terms  of  storage.  In  order  to  accurately  measure  the  efficiency 
of  these  methods,  we  define  the  effective  SSP  coefficient  of  a  method  as  Ce//  = 
this  iiormalizatioii  enables  us  to  compare  the  cost  of  integration  up  to  a  given  time 
using  diagonally  implicit  schemes  of  order  p  >  1.  Unfortunately,  the  optimal  implicit 
SSP  methods  have  effective  SSP  coefficient  less  than  or  equal  to  two,  making  them 
probably  too  inefficient  for  practical  use.  \\^  list  effective  SSP  coefficients  of  the 
numerically  optimal  methods  in  Table  2.1.  The  coefficients  of  the  most  efficient 
representations  of  SSP  implicit  Runge-Kutta  methods  are  available  online  [11]. 

The  SSP  condition  provides  a  guarantee  of  other  necessary  properties.  When 
considering  implicit  Runge-Kutta  methods,  it  is  important  to  determine  whether 
there  exists  a  unique  solution  of  the  stage  equations.  The  strong  stability  preserving 
tiiiiestep  restriction  turns  out  to  be  sufficient  for  this  as  well  [14,  Theorem  7.1]. 
Furthermore,  the  SSP  condition  serves  to  guarantee  that  the  errors  introduced  in  the 
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rn  /  p 

2 

Implicit  Methods  | 

3  4  5  6 

1 

2 

- 

- 

- 

- 

2 

2 

1.37 

- 

- 

- 

3 

2 

1.61 

0.68 

- 

- 

4 

2 

1.72 

1.11 

0.29 

5 

2 

1.78 

1.21 

0.64 

6 

2 

1.82 

1.30 

0.83 

0.030 

7 

2 

1.85 

1.31 

0.89 

0.038 

8 

2 

1.87 

1.33 

0.94 

0.28 

9 

2 

1.89 

1.34 

0.99 

0.63 

10 

2 

1.90 

1.36 

1.01 

0.81 

11 

2 

1.91 

1.38 

1.03 

0.80 

Table  2.1:  Effective  SSP  coefficients  of  best  known  implicit  methods.  A  dash  indicates 
that  SSP  methods  of  this  type  cannot  exist.  A  blank  space  indicates  that  no  SSP 
methods  of  this  type  were  found. 


solution  of  the  stage  equations  due  to  numerical  roundoff  and  (for  implicit  methods) 
errors  in  the  implicit  solve  are  not  unduly  amplified  [14,  Theorem  7.2]. 

In  summary,  we  have  gone  further  than  we  proposed  or  anticipated  possible  in 
the  study  of  SSP  implicit  Runge-Kutta  methods.  We  have  found  optimal  methods  of 
order  up  to  six  and  up  to  eleven  stages,  which  are  diagonally  implicit  and  which  have 
sparse  representations,  th\is  making  them  more  efficient  for  implementation.  Wc  also 
have  enough  information  to  conjecture  that  the  optimal  effective  SSP  coefficient  over 
this  class  of  methods  is  bounded  by  Ce//  <  2. 

Ill  addition  to  our  results,  the  methodology  we  adopted  in  this  search  also  led  to 
work  that  is  beyond  the  scope  of  this  grant.  David  Kctcheson  has  used  the  ideas 
and  techniques  developed  in  the  proeess  of  this  rcseareh  to  find  low  storage  optimal 
explicit  SSP  Runge-Kutta  methods  of  order  up  to  four  and  of  many  stages  [13]. 

Year  3:  The  next  step  in  our  research  involved  the  preliminary  testing  of  implicit 
and  explicit  SSP  methods  on  a  variety  of  problems.  We  carried  out  many  numerical 
experiments  which  showed  the  need  for,  and  benefit  of  SSP  methods. 

Using  a  nonlinear  example,  we  showed  that  even  when  the  spatial  discretization 
is  total  variation  diminishing  (TVD)  when  coupled  with  forward  Euler  integration, 
this  is  not  sufficient  to  guarantee  that  it  will  be  TVD  when  combined  with  a  higher 
order  timc-discretization.  We  considered  Burgers’  equation  with  a  sine  wave  initial 
condition  and  periodic  boundary  conditions.  The  solution  is  right-travelling  and  over 
time  steepens  into  a  shock.  We  discretize  using  a  first  order  conservative  upwind 
approximation  which  is  TVD  for  At  <  Ax  when  coupled  with  forward  Euler.  Using 
this  fact  we  can  conclude  that  if  we  integrate,  instead,  using  backward  Euler,  the 
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solution  will  bo  TVD  for  all  values  of  At.  However,  when  coupled  with  second-order 
A-stable  implicit  trapezoidal  rule  or  the  A-stable,  L-stable,  and  B-stable  implicit 
midpoint  rule,  this  is  not  TVD  for  At  >  2Ax. 

Using  a  iion-SSP  explieit  Runge-Kutta  with  a  second  order  TVD  fiux-differonciiig 
method  with  the  superbee  slope  limiter,  we  further  demonstrated  that  the  tiniestep 
restriction  associated  with  the  linear  SSP  property  does  not  suffice  to  give  reasonably 
good  behavior  in  the  nonlinear  ease. 

We  also  performed  experiments  of  SSP  methods  coupled  with  the  weighted  essen¬ 
tially  lion-oscillatory  method.  We  observe  advantages  to  the  use  of  SSP  methods  for 
WENO  methods  on  linear  and  nonlinear  problems.  The  time-step  at  which  the  total 
variation  begins  to  rise  by  more  than  10”^^  is  much  higher  for  the  SSP  methods  than 
for  the  eorrespoiidiiig  noii-SSP  methods.  We  observe  that  in  each  ca^e  the  tiniestep 
restriction  for  L2  linear  stability  is  larger  than  that  required  for  the  TVD  property, 
and  that  the  non-SSP  method  is  less  efficient  than  the  SSP  methods. 

For  SSP  Runge-Kutta  methods,  it  is  desirable  that  the  internal  stages  also  be 
strongly  stable.  This  means  requiring  not  only  that  <  ||u’^||,  but  also  that 

each  stage  for  i  =  1,  ...,m  satisfy  Since  the  SSP  argument  relies 

on  convexity,  which  is  satisfied  at  the  intermediate  stages  as  well,  SSP  Runge-Kutta 
methods  have  intermediate  stage  SSP  properties.  The  SSP  guarantee  of  provable 
stability  even  for  the  intermediate  stages  is  given  with  no  additional  cost.  This 
condition  is  frequently  necessary  in  the  approximate  solution  of  hyperbolic  PDEs. 
For  example,  in  the  iiumerieal  solution  of  the  Euler  equations  of  gas  dynamics,  it  is 
important  that  negative  pressure  or  density  values  be  avoided  even  in  the  intermediate 
stages.  Violations  of  these  bounds  are  more  than  theoretically  problematic,  as  they 
lead  to  non-physical  states  and  typically  to  failure  of  the  solution  algorithm.  We 
considered  the  Rieniann  problem  for  the  Euler  equations  with  fifth-order  WEXO 
used  for  the  spatial-discretization.  When  we  determined  the  largest  CFL  number  a 
for  whieli  the  density  and  pressure  values  remain  positive  at  all  Runge-Kutta  stages, 
we  find  that  we  see  that  the  SSP  methods  allow  a  more  efficient  time-step  than  the 
iion-SSP  methods. 

We  examined  the  class  of  spectral  deferred  eorreetioii  methods  methods,  and 
demonstrated  that  they  can  be  written  as  explicit  Runge-Kutta  methods.  Using  this 
fact,  we  can  immediately  establish  bounds  on  the  SSP  coefficient  of  spectral  deferred 
correction  methods  and  also  eoiielude  that  downwind  operators  will  be  required  in 
order  for  explicit  spectral  DC  methods  to  be  SSP  if  they  are  of  order  greater  than 
four.  Similarly,  implicit  spectral  DC  methods  cannot  be  SSP  without  dowiiwinding 
if  their  order  exceeds  six. 

Finally,  David  Ketchesoii  independently  studied  the  SSP  properties  of  the  Runge- 
Kutta  Chebyshev  methods.  Verwers  second  order  methods  all  have  negative  Butcher 
eoe?eients.  so  they  are  not  SSP  under  any  positive  tiniestep.  He  found  first  and  second 
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order  SSP  methods  up  to  10  stages  that  have  the  theoretically  optimal  time-step. 
These  are  promising  for  fully  explicit  integration  of  convection-di?usion  equations 
without  operator  splitting.  Unlike  IMEX,  exponential  di?erencing,  etc.,  they  apply 
the  same  integration  method  to  the  sti?  and  nori-sti?  parts) 


3  Transitions 


Guowei  Wei  (Michigan  State  University)  and  Shan  Zhao  have  implemented  our  SSP 
methods  in  their  matched  interface  and  boundary  method  to  obtain  high  order  schemes 
in  both  space  and  time  for  hyperbolic  equations.  They  report  that  “Your  SSP  meth¬ 
ods  work  great!” . 

Zhilin  Li  at  North  Carolina  State  University  requested  the  coefficients  of  the  sec¬ 
ond  order  tow-stage  implicit  SSP  scheme  to  use  these  with  free  boundary/moving 
interface  problems  for  which  stability  is  always  an  issue.  I  was  able  to  advise  him  on 
how  to  apply  this  most  efficiently. 

Marsha  Berger  (NYU)  and  Uri  Shumlak  (University  of  Washington)  requested  the 
SSP  review  paper.  Additionally,  Marsha  Berger  requested  that  I  recommend  specific 
SSP  methods  from  the  paper. 

Francis  X.  Giraldo  (Naval  Postgraduate  School  in  Monterey,  CA)  contacted  me 
asking  about  the  theoretical  limits  on  the  order  of  SSP  explicit  methods. 
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5  Appendix  I 


The  optimal  implicit  Ruuge-Kutta  methods  found  under  this  grant 

Shu-Osher  Coefficients  for  Second  Order  methods:  The  optimal  s-stage  second 
order  implicit  SSP  Runge-Kutta  method  has  SSP  coefficient  2s  and  Shu-Osher  form 
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- 1 
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1 
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2s 

2s 

1  ■ 

•  = 
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•.  0 
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2s  J 

Shu-Osher  Coefficients  for  Third  Order  methods:  The  optimal  s-stage  third 
order  implicit  SSP  Runge-Kutta  method  has  SSP  coefficient  s  —  1  +  \/ \  and 
Shu-Osher  form 


where 
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'0 

/'2I 
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'•  /^ll 
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1  0 
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{s  4-  1)(6’  -  1  4-  -  1) 

s{s  +  1  +  —  1) 


Shn-Osher  Coefficients  for  Fourth  Order  methods; 


1.  Optimal  3-stage,  4th-order  method: 


/ill  =  0.157330905682085 
/i22  =  0.047573123554705 
1133  =  0.157021682372699 
/i42  =  0.079106848361263 
A21  =  0.703541497995214 
A41  =  0.168078141811591 
A43  =  0.549902549377947 


/i2i  =  0.342491639470766 
//32  =  0.338136048168635 
/i4i  =  0.081822264233578 
/i43  =  0.267698531248384 
A32  =  0.694594303739345 
A42  =  0.162500172803529 
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2.  Optimal  4-stage,  4tli-order  method: 


/ill  -  0.119309657880174 
/i2i  =  0.226141632153728  /i22 
//.32  =  0.180764254304414  //ga 
/i43  -  0.212545672537219  /i44 
/isi  -  0.010888081702583  /i52 
/i54  =  0.181099440898861  A21 
A32  =  0.799340893504885  A43 
A51  =  0.048147179264990  A52 
A54  =  0.800823091149145 

3.  Optimal  5-stage,  4th-order  method: 

//ii  =  0.072154507748981  /i2i 
/t22  =  0.071232036614272  /i32 
/i33  =  0.063186062090477  /i43 
/i44  =  0.077017601068238  ^54 
H55  =  0.106426690493882  pes 
Ai52  =  0.007472809894781  /i62 
A21  =  1  A32 

A43  =  0.934991917505507  A54 
Acs  =  0.894472670673021  A52 
A62  =  0.105527329326976 

4.  Optimal  6-stage,  4th-order  method: 

//•ii  =  0.077219435861458  P21 
H22  =  0.063842903854499  ^32 
/i33  =  0.058359965096908  p4i 
/i43  =  0.103230521234296 
/i54  =  0.128204308556197  ^33 
fiQ3  =  0.008043763906343  /igs 
/i66  =  0.077016336936138  ^73 
P76  =  0.114400114184912  A21 
A32  =  1  A41 

A43  =  0.805203213502341  A54 
A63  =  0.062741759593964  Aes 
A73  =  0.107673404480272  A76 


0.070605579799433 

0.070606483961727 

0.119309875536981 

0.034154109552284 

1 

0.939878564212065 

0.151029729585865 


0.165562779595956 

0.130035287184462 

0.154799860761964 

0.158089969701175 

0.148091381629243 

0.017471397966712 

0.785413771753555 

0.954864191619538 

0.045135808380468 


0.128204308556198 

0.128204308556197 

0.008458154338733 

0.058105933032597 

0.064105484788524 

0.120160544649854 

0.013804194371285 

1 

0.065974025631326 

1 

0.937258240406037 

0.892326595519728 
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5.  Optimal  7-stage,  4th-order  method; 


/ill  =  0.081324471088377  /i2i 
^22  =  0.051065224656204  /i32 
//.33  -  0.036491713577701  //43 
/i44  =  0.037028821732794  ^54 
/i55  =  0.040474271914787  /igs 
/i66  =  0.061352000212100  /i73 
/i76  =  0.088170205242212  /i77 
/i83  -  0.001561606596621  /i87 
A21  =  1  A32 

A43  =  0.865661994183934  A54 
Aes  =  1  A73 

A76  =  0.810375930105481  A83 
A87  -  0.985647210475246 


0.108801609187400 

0.108801609187400 

0.094185417979586 

0.108801609187400 

0.108801609187400 

0.020631403945188 

0.080145231879588 

0.107240002590779 

1 

1 

0.189624069894518 

0.014352789524754 


6.  Optimal  8-stage,  4tli-order  method: 


/ill  =  0.080355939553359 
^22  =  0.054617345411549 
//33  =  0.039438131644116 
/i44  =  0.032427875074076 
/i54  =  0.083174746150582 
/i65  =  0.093742212796061 
/i76  =  0.093742212796061 
/i84  -  0.021977226754808 
/i88  =  0.055606577879005 
A21  =  1 
A43  =  1 

A54  =  0.887270992114641 
A76  =  1 

A87  =  0.765556774271797 


/i.2i  =  0.093742212796061 
/i32  =  0.093742212796061 
//43  =  0.093742212796061 
/i.51  =  0.004426522032754 
/i55  =  0.030116385482588 
/i66  =  0.038334326442344 
/i77  =  0.058861620081910 
/i87  =  0.071764986041253 
//98  =  0.093742212796061 

-^32  =  1 

A51  =  0.047220157287989 

-^65  =  1 

A84  =  0.234443225728203 
A98  =  1 
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7.  Optimal  9-stage,  4tli-order  method: 


/^ll  = 

0.068605696784244 

/^21 

1^22  = 

0.048685583036902 

/^32 

/^33  = 

0.039925150083662 

(H3 

/i44  = 

0.031928917146492 

/^54 

/^55  = 

0.029618614941264 

/^61 

/^G2  — 

0.001326570052113 

/^65 

/^GG  = 

0.029699905991308 

/^76 

/^77  = 

0.035642110881905 

/^87 

/^88  = 

0.050978240433952 

/^95 

/i98  = 

0.065270626421385 

/^99 

MiO,9  = 

=  0.083046524401968 

A21 

A32  = 

0.936520713898770 

^^43 

^^54  = 

1 

'^G2  = 

0.015973817828813 

'^GS 

A7G  = 

1 

00 

A95  = 

A  10,9  = 

0.214047464461523 

=  1 

(X 

8.  Optimal  10-stage,  4th-order  method: 


/7ii  =  0.053637857412307 

/^21  = 

/X22  =  0.042472343576273 

/^32  = 

/i33  =  0.039816143518898 

/^43  = 

/744  =  0.034233821696022 

/^54  — 

//,55  =  0.030626774272464 

/^G5  = 

/icG  =  0.029485772863308 

1^12  — 

/i7G  =  0.064406146499568 

/i77  = 

/i87  =  0.073302847899924 

/^88  = 

/798  =  0.073302847899924 

/^99  = 

//io,6  =  0.012892211367605 

/^10,9  = 

/iio!io  =  0.053275700719583 

A21  =  1 

A32  = 

A43  =  0.990280128291965 

II 

'^GS  =  1 

A72  = 

A70  =  0.878630890132646 

As?  = 

Ags  =  1 

Aio,6  = 

Aio.g  =  0.824124004224143 

^11,10 

0.082269487560004 

0.077774790319743 

0.083046524401968 

0.083046524401968 

0.008747971137402 

0.072971983212453 

0.083046524401968 

0.083046524401969 

0.017775897980583 

0.057552171403649 

0.990643355064403 

1 

0.105338196876962 

0.878687985294225 

1 

0.785952535538477 


0.073302847899924 
0.063734820131903 
0.072590353622503 
0.073302847899924 
0.073302847899924 
0.008896701400356 
0.033369849008191 
0.037227578299133 
0.046126339053885 
=  0.060410636532319 
=  0.073302847899924 
0.869472632481021 
1 

0.121369109867354 

1 

-  0.175875995775857 
=  1 
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9.  Optimal  11-stage,  4tli-order  method; 


/ill  =  0.056977945207836 

H21  =  0.065880156369595 

/i22  =  0.043484869703481 

/i32  =  0.065880156369595 

//,33  -  0.035790792116714 

//41  =  0.000026595081404 

/M3  -  0.061212831485396 

/244  =  0.029306212740362 

/X54  =  0.065880156369595 

p,55  =  0.028274789742965 

/265  =  0.065880156369595 

Hm  =  0.025442782369057 

/i76  =  0.065880156369595 

H77  =  0.029602951078198 

/f83  =  0.009935800759662 

Hs7  =  0.055944355609932 

Uss  =  0.027887296332663 

Hqs  =  0.065880156369595 

/igg  =  0.033340440672342 

Pio,g  =  0.065880156369595 

/iio,io  =  0.042024506703707 

Huj  =  0.012021727578515 

/in’io  =  0.053858428791080 

/ill’ll  =  0.045164424313434 

HUM  =  0.065880156369595 

A21  =  1 

-^32  =  1 

A41  =  0.000403688802047 

A43  =  0.929154313811668 

A54  =  1 

^65  =  1 

A70  =  1 

A83  =  0.150816289869158 

A87  =  0.849183710130842 

Ag8  =  1 

Aio.o  =  1 

All, 7  =  0.182478734735714 

All, 10  =  0.817521265264286 

■^12,11  =  1 

Oshcr  Coefficients  for  Fifth  Order 

methods: 

Optimal  4-stage,  5tli-order  method: 

H2i  =  0.125534208080981 

M22 

=  0.125534208080983 

H32  =  0.350653119567098 

M33 

=  0.048181647388277 

(Ml  =  0.097766579224131 

/^42 

=  0.000000005345013 

=  0.404181556145118 

M44 

=  0.133639210602434 

H51  =  0.022869941925234 

M52 

=  0.138100556728488 

M53  =  0.157510964003014 

M54 

=  0.277310825799681 

A21  =  0.143502249669229 

A32 

=  0.400843023432714 

A41  =  0.111760167014216 

^^42 

=  0.000000006110058 

A43  =  0.462033126016285 

=  0.026143376902960 

A52  =  0.157867252871240 

^^53 

=  0.180055922824003 

A54  =  0.317003054133379 
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2.  Optimal  5-stage,  5th-order  method; 


=  0.107733237609082 
fi{2, 2)  =  0.107733237609079  /j(3,  1) 
//(3, 2)  =  0.205965878618791  //(3,  3) 
/i(4, 1)  =  0.010993335656900  ^^(4,  2) 
/i(4, 3)  =  0.245761367350216  ^(4, 4) 
/i(5, 1)  -  0.040294985548405  ^(5,  2) 
//(5, 3)  =  0.024232322953809  m(5,  4) 
/j(5, 5)  =  0.098999612937858  /i(6, 3) 
/v.(6, 4)  =  0.023678103998428  ^(6, 5) 
A(2, 1)  =  0.344663606249694  A(3, 1) 
A(3, 2)  =  0.658932601159987  A(4, 1) 
A(4, 2)  =  0.000000100208717  A(4, 3) 
A(5, 1)  =  0.128913001605754  A(5,2) 
A(5, 3)  -  0.077524819660326  A(5, 4) 
A(6, 3)  =  0.255260385110718  A(6, 4) 
A(6,5)  =  0.623567413728619 

3.  Optimal  6-stage,  5tli-order  method: 

fi{2, 1)  =  0.084842972180459  n{2,  2) 
//(3, 2)  =  0. 149945333907731  ^(3, 3) 
//(4, 3)  =  0.175767531234932  fi{4, 4) 
p(5, 1)  =  0.024709139041008  ^(5, 4) 
//(5, 5)  =  0.054767418942828  /7.(6,  2) 
//(6, 3)  =  0.026804592504486  /y(6,  5) 
/i(6, 6)  =  0.085074359110886  ^{7, 3) 
//(7, 4)  =  0.042600565019890  /i(7, 6) 
A(2, 1)  =  0.422021261021445  A(3, 2) 
A(4, 3)  =  0.874293218071360  A(5, 1) 
A(5, 4)  =  0.861728690085026  A(6,  2) 
A(6, 3)  =  0.133329934574294  A(6,  5) 
A(7, 3)  =  0.024117294382203  A(7, 4) 
A(7,6) =  0.752865185365536 


0.000009733684024 

0.041505157180052 

0.000000031322743 

0.079032059834967 

0.011356303341111 

0.220980752503271 

0.079788022937926 

0.194911604040485 

0.000031140312055 

0.035170229692428 

0.786247596634378 

0.036331447472278 

0.706968664080396 

0.075751744720289 


0.084842972180464 

0.063973483119994 

0.055745328618053 

0.173241563951140 

0.014574431645716 

0.159145416202648 

0.004848530454093 

0.151355691945479 

0.745849859731775 

0.122906844831659 

0.072495338903420 

0.791612404723054 

0.211901395105308 
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4.  Optimal  7-stage,  5tli-order  method: 


=  0.077756487471956  P22 
P32  =  0.126469010941083  //as 
//43  =  0.143639250502198  //44 
Psi  -  0.011999093244164  ^54 
P55  =  0.047108760907057  pg2 
Pg3  =  0.027138257330487  /xgs 
Pgg  =  0.037306165750735  /i73 
/i76  =  0.140855998083160  /i77 
Hni  =  0.009653207936821  pss 
PsG  =  0.000177781270869  ps? 
A21  =  0.482857811904546  A32 
A43  =  0.891981318293413  A51 
A54  =  0.900717090387559  Ag2 
Ag3  =  0.168525096484428  Ags 
A73  =  0.125302322168346  A7G 
A84  =  0.059945182887979  Ags 
Asg  =  0.001103998884730  A87 

5.  Optimal  8-stage,  5tli-order  method: 

/i2i  =  0.068228425119547  P22 
/i32  =  0.105785458668142  1133 
fi4i  =  0.119135238085849  (144 
P51  ^  0.009164078944895  H54 
//55  =  0.039406904101415  //g2 
/iG3  =  0.019703233696280  //g5 
Pgg  =  0.045239659320409  ^73 
Hn  =  0.116977452926909  (177 
/i84  =  0.011255581082016  /igs 
//87  =  0.114515518273119  //gs 
/i95  =  0.002607774587593  poG 
p.j8  =  0.104666894951906  A21 
A32  =  0.799508082567950  A43 
A51  -  0.069260513476804  A54 
Ag2  =  0.056144626483417  Ag3 
Agg  =  0.794939486396848  A73 
A7G  =  0.884095226988328  A84 
Ass  =  0.049438833770315  A87 
A95  =  0.019709106398420  A9G 
A98  =  0.791054172708715 


0.077756487471823 

0.058945597921853 

0.044443238891736 

0.145046006148787 

0.011454172434127 

0.122441492758580 

0.020177924440034 

0.077972159279168 

0.025430639631870 

0.124996366168017 

0.785356333370487 

0.074512829695468 

0.071128941372444 

0.760345962143127 

0.874697677831654 

0.157921009644458 

0.776211398253764 


0.068228425081188 

0.049168429086829 

0.040919294063196 

0.120257079939301 

0.007428674198294 

0.105180973170163 

0.015335646668415 

0.050447703819928 

0.006541409424671 

0.060382824328534 

0.024666705635997 

0.515658560550227 

0.900403391614526 

0.908882077064212 

0.148913610539984 

0.115904148048060 

0.085067722561958 

0.865488353423280 

0.186426667470161 
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6.  Optimal  9-stage,  5tli-order  method: 


/i2i  =  0.057541273792734 
H-n  =  0.089687860942851 
//,43  =  0.101622955619526 
//5i  =  0.009276188714858 
//55  =  0.040815264589441 
//65  -  0.101125244372555 
Hn  =  0.003606182878823 
P76  =  0.090586614534056 
=  0.011070977346914 
//88  =  0.046669302312152 
=  0.102117191974435 
/xio,6  =  0.000157554758807 
/xio’g  =  0.088454624345414 
A32  =  0.797947256574797 
A51  =  0.082529667434119 
A(i2  =  0.100295062538531 
A73  =  0.032083982209117 
A76  =  0.805943410735452 
A87  =  0.901502211016037 
A98  =  0.908530232837680 
Aio,7  =  0.210035759124536 


^22  =  0.057541282875429 
//33  -  0.041684970395150 
//44  =  0.040743690263377 
//54  =  0.101958242208571 
/i62  =  0.011272987717036 
//66  =  0.040395338505384 
im  =  0.018205434656765 
Ai77  =  0.042925976445877 
//87  =  0.101327254746568 
Ai95  =  0.010281040119047 
Hgg  =  0.050500143250113 
//io,7  =  0.023607648002010 
A21  =  0.511941093031398 
A43  =  0.904133043080300 
A54  =  0.907116066770269 
Acs  =  0.899704937426848 
A74  =  0.161972606843345 
A84  =  0.098497788983963 
A95  =  0.091469767162319 
Aio,6  =  0.001401754777391 
Aio’g  =  0.786975228149903 
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7.  Optimal  10-stage,  5th-order  method: 


//2i  =  0.052445615058994 
/i32  =  0.079936220395519 
//43  =  0.089893189589075 
/isi  =  0.007606429497294 
/i55  =  0.035536573874530 
/X65  =  0.089447242753894 
//73  =  0.003271387942850 
P76  =  0.080215515252923 
//84  -  0.009638972523544 
m  =  0.040785658461768 
p.98  =  0.089540979697808 
/iio.c  =  0.005634796609556 
/iio,9  =  0.086547180546464 
//.117  =  0.001872759401284 
/til, 10  =  0.079160150775900 
^2  =  0.809542670828687 
A51  =  0.077033029836054 
A62  =  0.094135396158718 
A73  =  0.033130514796271 
A76  =  0.812371189661489 
A87  =  0.902382678155958 
A98  =  0.906813500744962 
Aio,7  =  0.066440169285130 
All, 7  =  0.018966103726616 
All, 10  =  0.801683136446066 


H22  =  0.052445635165954 
//33  =  0.038724845476313 
im  =  0.037676214671832 
/i54  =  0.090180506502554 
P<32  =  0.009295158915663 
//6G  =  0.036490114423762 
P74  =  0.015255382390056 
//77  =  0.035768398609662 
1187  =  0.089103469454345 
/i95  =  0.009201462517982 
yigg  =  0.042414168555682 
/iio,7  =  0.006560464576444 
Atio.io  =  0.043749770437420 
//i’i,8  =  0.017616881402665 
A21  =  0.531135486241871 
A43  =  0.910380456183399 
A54  =  0.913290217244921 
Acs  =  0.905864193215084 
A74  =  0.154496709294644 
A84  =  0.097617319434729 
A95  =  0.093186499255038 
Aio,g  =  0.057065598977612 
Aio’9  =  0.876494226842443 
All, 8  =  0.178412453726484 
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8.  Optiirial  11-stage,  5th-order  method; 


/X21  =  0.048856948431570 
/X32  -  0.072383163641108 
//.43  =  0.080721632683704 
/X51  =  0.006438090160799 
/X55  =  0.032672027896742 
//63  =  0.000719846382100 
/X66  =  0.033437798720082 
P74  =  0.012192534706212 
P77  =  0.033377699686911 
/X87  =  0.079986775597087 
^95  =  0.008095394925904 
p.99  =  0.036372965664654 
/xio,7  =  0.005394911565057 
/iio.io  =  0.032282094274356 
/Xu, 8  =  0.008920593887617 
/xii,n  =  0.042478561828713 
/xi2,9  =  0.011637432775226 
A21  =  0.553696439876870 
A43  =  0.914819326070196 
A54  =  0.918370981510030 
Ag3  =  0.008158028526592 
A73  =  0.034327672500586 
A76  =  0.827494171134198 
A87  =  0.906491181031666 
Agg  -  0.908254782302260 
Aio,7  =  0.061140603801867 
Aii’7  =  0.040471104837131 
An,io  =  0.858431687176596 
Ai2,9  =  0.131887178872293 


/X22  =  0.048856861697775 
/X33  =  0.035920513887793 
/X44  =  0.034009594943671 
/X54  =  0.081035022899306 
/X62  =  0.007591099341932 
/X65  =  0.079926841108108 
/X73  =  0.003028997848550 
/X76  -  0.073016254277378 
/x84  =  0.008251011235053 
/X88  =  0.035640440183022 
/xgs  =  0.080142391870059 
/xio.e  =  0.005907318148947 
=  0.076935557118137 
/xiij  =  0.003571080721480 
/iii,io  =  0.075746112223043 
/xi2,8  =  0.004170617993886 
/xi2,n  =  0.072377330912325 
^2  =  0.820319346617409 
A51  =  0.072962960562995 
A62  =  0.086030028794504 
Aes  =  0.905811942678904 
A74  =  0.138178156365216 
A84  =  0.093508818968334 
Ays  =  0.091745217287743 
Aio,6  =  0.066947714363965 
Aio’g  =  0.871911681834169 
An^  =  0.101097207986272 
A123  =  0.047265668639449 
Ai2,ii  =  0.820253244225314 
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Shu-Osher  Coefficients  for  sixth  Order  methods: 


1.  Tlic  optimal  sixth-order,  six-stage  method  (c  =  0.18): 


//2i  =  0.306709397198437 
/i3i  =  0.100402778173265 
/X33  =  0.100402700098726 
/i42  =  0.000708584139276 
//.44  -  0.028228318307509 
1152  =  0.000026687930165 
H54  =  0.331296656179688 
Hai  =  0.000033015066992 
/i63  =  0.395057247524893 
^65  =  0.421912313467517 
liji  =  0.054129307323559 
=  0.233976271277479 
/X75  =  0.303060566272042 
A21  =  0.055928810359256 
A32  =  0.000000002666388 
A42  =  0.000129211130507 
A51  =  0.018587746937629 
A53  =  0.024929494718837 
Aei  =  0.000006020335333 
Asa  =  0.072039142196788 
Ass  =  0.076936194272824 
A72  =  0.000379944400556 
A74  =  0.033716209818106 
A76  =  0.024795346049276 


H22  =  0.306709397198281 
^732  =  0.000000014622272 
/X41  =  0.000015431349319 
^l43  =  0.383195003696784 
fi5i  =  0.101933808745384 
fi53  =  0.136711477475771 
fi55  =  0.107322255666019 
P62  =  0.000000017576816 
fiM  =  0.014536993458566 
PsG  =  0.049194928995335 
fi72  =  0.002083586568620 
P74  =  0.184897163424393 
fire  =  0.135975816243004 
A31  =  0.018308561756789 
A41  =  0.000002813924247 
A43  =  0.069876048429340 
A52  =  0.000004866574675 
A54  =  0.060412325234826 
A62  =  0.000000003205153 
A64  =  0.002650837430364 
Ati  =  0.009870541274021 
A73  =  0.042665841426363 
A75  =  0.055263441854804 
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2.  The  optimal  sixth-order,  seven-stage  method  (c  =  0.26): 


//2i  -  0.090485932570398 
^22  =  0.090485932570397 
/i32  -  0.346199513509666 
/i33  =  0.056955495796615 
//4i  -  0.089183260058590 
/i42  =  0.122181527536711 
/r43  =  0.340520235772773 
1^44  =  0.086699362107543 
/isi  =  0.214371998459638 
fi52  =  0.046209156887254 
/i53  =  0.215162143673919 
(154  =  0.000000362542364 
(155  =  0.209813410800754 
p.6i  =  0.000000591802702 
(152  =  0.390556634551239 
P63  =  0.000000491944026 
(154  =  0.330590135449081 
(155  =  0.007410530577593 
(155  =  0.070407008959133 
(171  =  0.000000021842570 
(172  =  0.325421794191472 
(173  =  0.069025907032937 
(1.74  =  0.373360315300742 
(175  =  0.007542750523234 
(175  =  0.005465714557738 
(177  =  0.063240270982556 
(isi  =  0.044161355044152 
(182  =  0.204837996136028 
(i83  =  0.191269829083813 
(184  =  0.255834644704751 
(i85  =  0.015984178241749 
(185  =  0.016124165979879 
(187  =  0.151145768228502 


A21  =  0.023787133610744 
A32  =  0.091009661390427 
A41  =  0.023444684301672 
A42  =  0.032119338749362 
A43  =  0.089516680829776 
A51  =  0.056354565012571 
A52  =  0.012147561037311 
A53  =  0.056562280060094 
A54  =  0.000000095305905 
Aei  =  0.000000155574348 
A62  =  0.102670355321862 
A63  =  0.000000129323288 
A64  =  0.086906235023916 
Aes  -  0.001948095974350 
A71  =  0.000000005742021 
A72  =  0.085547570527144 
A73  =  0.018145676643359 
A74  =  0.098149750494075 
A75  =  0.001982854233713 
A76  =  0.001436838619770 
Agi  =  0.011609230551384 
A82  =  0.053848246287940 
A83  =  0.050281417794762 
A84  =  0.067254353278777 
A85  =  0.004201954631994 
A86  =  0.004238754905099 
A87  =  0.039733519691061 
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3.  The  optimal  sixth-order,  eight-stage  method  (c  =  2.25); 


/r2i  =  0.078064586430339 
fi22  =  0.078064586430334 
P31  =  0.000000000128683 
/i32  =  0.207887720440412 
/i33  =  0.051491724905522 
mi  =  0.039407945831803 
P43  =  0.256652317630585 
1144  =  0.062490509654886 
/i5i  -  0.009678931461971 
m2  =  0.113739188386853 
/i54  =  0.227795405648863 
//55  =  0.076375614721986 
m2  =  0.010220279377975 
m3  =  0.135083590682973 
//(i5  =  0.235156310567507 
ma  =  0.033370798931382 
m2  =  0.000000009428737 
//73  =  0.112827524882246 
P74  =  0.001997541632150 
/i75  =  0.177750742549303 
ma  =  0.099344022703332 
mi  =  0.025183595544641 
mi  =  0.122181071065616 
m2  =  0.000859535946343 
ms  =  0.008253954430873 
ms  =  0.230190271515289 
ms  =  0.046429529676480 
P86  =  0.017457063072040 
mi  =  0.017932893410781 
ms  =  0.322331010725841 
mi  =  0.011069087473717 
m2  =  0.010971589676607 
ms  =  0.068827453812950 
ms  =  0.048864283062331 
ms  =  0.137398274895655 
ms  =  0.090347431612516 
mi  =  0.029504401738350 
ms  =  0.000167109498102 


A21  =  0.175964293749273 
A31  =  0.000000000290062 
A32  =  0.468596806556916 
A41  =  0.088828900190110 
A43  =  0.578516403866171 
A51  =  0.021817144198582 
A52  =  0.256377915663045 
A54  =  0.513470441684846 
A62  =  0.023037388973687 
A63  =  0.304490034708070 
Ao5  =  0.530062554633790 
A72  =  0.000000021253185 
A73  =  0.254322947692795 
A74  =  0.004502630688369 
A75  =  0.400665465691124 
A7G  =  0.223929973789109 
Agi  =  0.275406645480353 
A82  =  0.001937467969363 
A83  =  0.018605123379003 
A84  =  0.518868675379274 
Ags  =  0.104656154246370 
Age  =  0.039349722004217 
A87  =  0.040422284523661 
Agi  =  0.024950675444873 
A92  =  0.024730907022402 
A93  =  0.155143002154553 
A94  =  0.110144297841125 
A95  =  0.309707532056893 
Age  =  0.203650883489192 
A97  =  0.066505459796630 
Agg  =  0.000376679185235 
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4.  The  optimal  sixth-order,  nine-stage  method  (c  =  5.80): 


P21 

P22 

M31 

/^32 

/-i33 

/^42 

M43 

P44 

/'51 

At52 

/t53 

/i54 

M55 

f>62 

Pgs 

/'G4 

/^G5 

P71 

/'72 

A'73 

/i75 

/^7G 

/'77 

/'81 

/^82 

/^83 

/-f84 

//85 

P86 

M87 

/'88 

M'Jl 

^92 

P93 

/'94 

P95 

/^9G 

P97 

M98 

M99 

MlO,l 

/^10,2 

M10,3 

/ilO.4 

^10,5 


0.060383920365295 

0.060383920365140 

0.000000016362287 

0.119393671070984 

0.047601859039825 

0.000000124502898 

0.144150297305350 

0.016490678866732 

0.014942049029658 

0.033143125204828 

0.020040368468312 

0.095855615754989 

0.053193337903908 

0.000006536159050 

0.000805531139166 

0.015191136635430 

0.054834245267704 

0.089706774214904 

0.000006097150226 

0.018675155382709 

0.025989306353490 

0.000224116890218 

0.000125522781582 

0.125570620920810 

0.019840674620006 

0.000000149127775 

0.000000015972341 

0.034242827620807 

0.017165973521939 

0.000000000381532 

0.001237807078917 

0.119875131948576 

0.056749019092783 

0.000000072610411 

0.000000387168511 

0.000400376164405 

0.000109472445726 

0.012817181286633 

0.011531979169562 

0.000028859233948 

0.143963789161172 

0.060174596046625 

0.001577092080021 

0.000008909587678 

0.000003226074427 

0.000000062166§f0 

0.009112668630420 

0.008694079174358 


A21 

-^31 

-^32 

-^42 

-^43 

A5I 

A52 

A53 

A54 

Aei 

Ag2 

^63 
^64 
A  65 

A7I 
A72 
A73 
A74 
A  75 
A76 
Asi 

A  82 

A83 
A84 
A  85 

A86 

A87 

A9I 

A92 

A93 


A94 

A95 

A9G 

A97 

A98 

AlO.l 

Ai0,2 

Ai0,3 

Ai0,4 

Ai0,5 

Ai0,6 

Aloj 

Ai0,8 

A  10,9 


0.350007201986739 

0.000000094841777 

0.692049215977999 

0.000000721664155 

0.835547641163090 

0.086609559981880 

0.192109628653810 

0.116161276908552 

0.555614071795216 

0.000037885959162 

0.004669151960107 

0.088053302494510 

0.317839263219390 

0.519973146034093 

0.000035341304071 

0.108248004479122 

0.150643488255346 

0.001299063147749 

0.000727575773504 

0.727853067743022 

0.000000864398917 

0.000000092581509 

0.198483904509141 

0.099500236576982 

0.000000002211499 

0.007174780797111 

0.694839938634174 

0.000000420876394 

0.000002244169749 

0.002320726117116 

0.000634542179300 

0.074293052394615 

0.066843552689032 

0.000167278634186 

0.834466572009306 

0.009141400274516 

0.000051643216195 

0.000018699502726 

0.000000360342058 

0.052820347381733 

0.050394050390558 

0.103597678603687 

0.159007699664781 

0.624187175011814 


5.  The  optimal  sixth-order,  ten-stage  method  (c  =  8.10); 


//2i  =  0.054638144097621 
;222  =  0.054638144097609 
/i32  =  0.094708145223810 
/i33  =  0.044846931722606 
/i43  =  0.108958403164940 
=  0.031071352647397 
^51  =  0.004498251069701 
/i52  -  0.005530448043688 
/i54  =  0.107851443619437 
/i55  =  0.018486380725450 
/i(i2  =  0.015328210231111 
Pgs  =  0.014873940010974 
P64  =  0.000000013999299 
P65  =  0.093285690103096 
/i66  =  0.031019852663844 
/i73  =  0.023345108682580 
P74  =  0.000000462051194 
//76  =  0.100142283610706 
P77  =  0.037191650574052 
/284  =  0.020931607249912 
/i85  =  0.007491225374492 
/i86  =  0.000000004705702 
//.87  =  0.094887152674486 
/i88  =  0.041052752299292 
/ig4  -  0.000000000437894 
/i95  =  0.013484714992727 
/ige  =  0.012301077330264 
/igg  -  0.097178530400423 
//.gg  =  0.039273658398104 
/iio,i  =  0.000987065715240 
/iiog  =  0.000000347467847 
/iio'c  =  0.004337021151393 
/iioj  =  0.011460261685365 
//.103  =  0.002121689510807 
/xio’g  =  0.104338127248348 
/iio.io  =  0.042268075457472 
/ill, 3  =  0.000656941338471 
Miij  =  0.015039465910057 
Mips  =  0.004816543620956 
//ai,g  -  0.031302441038151 
10  =  0.071672462436845 


A21  -  0.442457635916190 
A32  =  0.766942997969774 
A43  =  0.882341050812911 
A51  =  0.036426667979449 
A52  =  0.044785360253007 
A54  =  0.873376934047102 
A62  =  0.124127269944714 
A63  =  0.120448606787528 
A64  =  0.000000113365798 
Acs  =  0.755424009901960 
A73  =  0.189047812082446 
A74  =  0.000003741673193 
A76  =  0.810948446244362 
A84  =  0.169503368254511 
Ass  =  0.060663661331375 
Ase  =  0.000000038106595 
A87  =  0.768392593572726 
Ag4  =  0.000000003546047 
A95  =  0.109198714839684 
Agg  =  0.099613661566658 
Ags  =  0.786948084216732 
Aio,i  =  0.007993221037648 
Aio’2  =  0.000002813781560 
Aio’o  =  0.035121034164983 
Aioj  =  0.092804768098049 
Aio, 8  =  0.017181361859997 
Aio’g  =  0.844926230212794 
All’s  =  0.005319886250823 
All, 7  =  0.121789029292733 
A113  =  0.039004189088262 
Aii’g  =  0.253485990215933 
All, 10  =  0.580400905152248 
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6  Appendix  2 


This  appendix  contains  MATLAB  scripts  which  we  created  as  part  of  this  project. 

•  Evaluate  the  Radius  of  absolute  rnonotonicity  from  the  Butcher  array: 

function  r  =  am_radius(A,b) 

/oBy  David  Ketcheson 

"/oEvaluates  the  Radius  of  absolute  rnonotonicity 
7oOf  a  Runge-Kutta  method,  given  the  Butcher  array. 

1 

7oFor  an  m-stage  method,  A  should  be  cui  m  x  m  matrix 
7oand  b  should  be  a  column  vector  of  length  m. 

7o 

7oAccuracy  can  be  changed  by  modifying  the  value  of  eps. 

7oMethods  with  very  large  radii  of  a.m.  (>50)  will  require 
7,rmax  to  be  increased. 

rmax=50 ;  eps=l . e-12 ; 

m=length(b) ;  e=ones(m, 1) ; 

K=  [A;b^] ; 
rlo=0;  rhi=rmax; 

while  rhi-rlo>eps  7oUse  bisection 
r=0.5*(rhi+rlo) ; 

X=eye(m)+r*A;  beta=K/X;  ech=r*K*(X\e) ; 
if  (min(beta( : ))<-3.e-16  | I  max(ech( : ) )>1 . +3 . e-16) 
rhi=r ; 
else 
rlo=r ; 
end 
end 

if  rhi==rmax  7o  r>=rmax 

error ( ’Error ;  increase  value  of  rmax  in  am_radius .m’ ) ; 
else 
r=rlo ; 
end 

•  Butcher  array  from  the  Shu  Osher  array: 
function  [A, b, c]=shuosher2butcher (alpha, beta) ; 
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7ofunction  [A,b,c]=shuosher2butcher(laiiibda,inu)  ; 

I 

7oBy  David  Ketcheson 
7o 

7oGenerate  Butcher  form  of  a  Runge-Kutta  method, 

7ogiven  its  Shu-Osher  or  modified  Shu-Osher  form 

y« 

yoFor  an  m-stage  method,  alpha  and  beta  (or  lambda  and  mu)  should  be 
7omatrices  of  dimension  (m+1)  x  m 
7o 

7oNote  that  MATLAB  indexes  from  1,  while  the  Shu-Osher  coefficients 
7oare  usually  indexed  from  zero. 

s=size(alpha,2) ; 

X=eye(s)-alpha(l : end-1 ,  : )  ; 

A=X\beta(l : end-1 , : ) ; 

b=beta(end, : )+alpha(end, : )*A;  b=b^ ; 

c=sum(A,2) ; 

•  Butcher  to  modified  Shu-Osher:  This  function  generates  the  modified  Sliu-Osher 
form  of  a  Runge-Kutta  method,  given  its  Butcher  form  and  radius  of  absolute 
moriotonicity 

function  [lambda, mul] =butcher2modshuosher (A,b,r) ; 

7oBy  David  Ketcheson 

7oGenerate  modified  Shu-Osher  form  of  a  Runge-Kutta  method, 

7ogiven  its  Butcher  form  and  radius  of  absolute  monotonicity 
7o 

yoFor  an  m-stage  method,  A  should  be  an  m  x  m  matrix 
7oand  b  should  be  a  column  vector  of  length  m. 

7o 

7oNote  that  MATLAB  indexes  from  1,  while  the  Shu-Osher  coefficients 
7oare  usually  indexed  from  zero. 

Aup=triu(A) ; 

if  max(abs(Aup))>0  mclass= ' implicit ' ;  else  mclass= ' explicit ' ;  end 
if  nargin<3  r  =  am_radius(A,b) ;  end 

s=size(A, 1) ; 

K=  [A;b’] ; 

G=eye(s)+r*A; 
mul=K/G; 
lambda=r*mul ; 
for  i=l:s+l 

if  strcmpCmclass,  ^implicit O  7oO  stage  is  u_n 
alpha(i)=l-sum(lambda(i , 1 : s) ) ; 
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end 

end 

"/oEliminate  diagonal  terms  of  lambda  array  (for  implicit  methods) 
lambda=lambda;  mul=mul ; 
for  i=l:s 

if  larabdaCi  ,  i)  "'=1 

f ac=l/( 1-lambda (i,i)) ; 
mul(i, :)=fac*mul(i, :) ; 
lambdaCi , : )=f ac*lambda(i , : ) ; 
lambdaCi , i)=0 . ; 
end 
end 


•  Butcher  to  Shu-Osher: 

function  [alpha,beta]=butcher2shuosher(A,b,r) ; 

7oBy  David  Ketcheson 
7o 

7oGenerate  Shu-Osher  form  of  cin  explicit  Runge-Kutta  method, 

7ogiven  its  Butcher  form  and  radius  of  absolute  monotonicity 
7o 

7oFor  an  m-stage  method,  A  should  be  cin  m  x  m  matrix 
7oand  b  should  be  a  column  vector  of  length  m. 

7o 

7oNote  that  MATLAB  indexes  from  1,  while  the  Shu-Osher  coefficients 
7oare  usually  indexed  from  zero. 

if  narginO  r  =  am_radius(A,b) ;  end 

s=size(A, 1) ; 

K=[A;b']  ; 

G=eye(s)+r*A; 
beta=K/G; 
alpha=r*beta; 
for  i=2:s+l 

alpha (i , l)=l“sum(alpha(i ,2 : s)) ; 
end 

•  Radius  of  circle  contractivity  of  a  Runge-Kutta  method^  given  the  Butcher  army: 


function  [r]  =  cc_radius(A,b) 
7oBy  David  Ketcheson 
7o 
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VoEvaluates  the  Radius  of  circle  contractivity 
"/oOf  a  Runge-Kutta  method,  given  the  Butcher  array 
7o 

7oFor  an  m-stage  method,  A  should  be  an  m  x  m  matrix 
7oand  b  should  be  a  column  vector  of  length  m. 

7o 

"/oAccuracy  can  be  changed  by  modifying  the  value  of  eps . 

VoMethods  with  very  large  radii  of  a,m,  (>1000)  will  require 
7ormax  to  be  increased, 

rmax=1000;  eps=l , e-13; 

if  min(b)<=0 
r=0; 
else 

m=length(b) ; 

B=diag(b) ; 

M=B*A+A'*B-b*b' ; 
rlo=0;  rhi=rmax; 
while  rhi“rlo>eps 
r=0 . 5*(rhi-i-rlo) ; 

X=M+B/r ; 

if  min(eig(X))<“3.e“16 
rhi=r ; 
else 
rlo=r ; 
end 
end 
end 

if  rhi==rmax  7  r>=rmax 

error ( ’Error :  increase  value  of  rmax  in  cc_radius .m’ ) ; 
else 
r=rlo; 
end 

•  Make  the  Butcher  arrays  for  different  methods: 

function  [A , b , c , r] =makebutcher (name , s) 

"/oBy  David  Ketcheson 

7 

7Set  up  Butcher  arrays  A,b,c  for  various  methods 
7A1so  returns  SSP  coefficient  r 

7For  families  of  methods,  optional  input  s  is  the  number  of  stages 
if  nargin<2  s=l;  end 
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switch  name 


*/.=================SSP  Methods=======^ 

*/,=================Explicit  Methods 


case  'FEll’ 

"/.Forward  Euler 
s=l;  r=l; 

A=[0]  ; 

b=[l]  '  ;  c=[0]  ’  ; 

case  'SSP22' 
s=2;  r=l; 

A=[0  0;  10]; 

b=[l/2  1/2]';  c=sum(A,2); 

case  'SSP42’ 
s=4;  r=3 ; 

A=[0  0  0  0;  1/3000;  1/31/300;  1/3  1/3  1/3  0] ; 
b=l/4*ones (m, 1) ;  c=sum(A,2); 

case  'SSP33' 
s=3;  r=l; 

A=[0  0  0;  100;  1/4  1/4  0] ; 
b=[l/6  1/6  2/3]';  c=sum(A,2); 

case  'SSP43' 
s=4;  r=2; 

A=[0  000;  1/200  0;  1/2  1/2  0  0;  1/6  1/6  1/6  0]; 
b=[l/6  1/6  1/6  1/2]';  c=sum(A.2); 

case  ’SSP104' 
s=10;  r=6; 

alphaO=diag(ones (1 , s-1) , -1) ; 
alphaO (6 . 5) =2/5 ;  alphaO (6 . 1) =3/5 ; 
betaO  =l/6*diag(ones(l ,s-l) ,-l) ; 
beta0(6.5)=l/15; 

A= (eye (s) -alphaO) \betaO ; 
b=l/10*ones(s, 1) ;  c=sum(A,2); 

case  ’rSSPs2' 

*/, Rational  (optimal,  low-storage)  s-stage  2nd  order  SSP 

if  s<2  error ( 'Explicit  second  order  SSP  family  requires  s>=2');  end 

r=s-l ; 

alpha= [zeros (1 , s) ; eye(s) ; ] ; 
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alpha(s+l, s)=(s-l)/s; 
beta=alpha/r ; 
alpha(s+l, l)=l/s; 

A=(eye(s) -alphaCl : s , : ) )\beta(l : s , : ) ; 
b=beta(s+l, : )+alpha(s+l , :)*A;  b=b^ ; 
c=siim(A,2) ; 

case  ^rSSPsS^ 

7oRational  (optimal,  low'-storage)  s"2“Stage  3rd  order  SSP 
if  roundCsqrt  (s) )  "'=sqrt  (s)  |  |  s<4 

error ( ^Explicit  third  order  SSP  family  requires  s=n^2,  n>lO; 

end 

n=s'‘2;  r=n“'s; 

alpha= [zeros ( 1 , n) ; eye (n) ; ] ; 

alphaCs* (s+l)/2+l,s*(s+l)/2)=(s-l)/(2*s-l) ; 

beta=alpha/r ; 

alpha(s*(s+l)/2+l, (s-l)*(s-2)/2+l)=s/(2*s-l) ; 

A=(eye(n) -alphaCl : n, : ) ) \beta(l :n, : ) ; 
b=beta(n+l , : )+alpha(n+l , : ) *A;  b=b^ ; 
c=sum(A,2) ; 


7,=================Implicit  Methods========================= 

case  'BEll' 

“/oBackward  Euler 
s=l;  r=l.elO; 

A=[l]; 

b=[l]^;  c=[l]^ 

case  ^SDIRK34'  7p3-stage,  4th  order  singly  diagonally  implicit  (SSP) 
s=3;  r=1.7588; 

g=0 . 5* (l-cos(pi/18)/sqrt (3) -sin (pi/ 18) ) ; 
q=(0.5-g)^2; 

A= [g  0  0 

0.5-g  g  0 

2*g  l’4*g  g] ; 

b=[l/(24*q)  l-l/(12*q)  l/(24*q)] ^ 
c=sum(A, 2) ; 

case  'ISSPm2^ 

7oOptimal  DIRK  SSP  schemes  of  order  2 
r=2*s; 

i=repmat ( (1 : s)  \  1 ,s) ;  j=repmat (1 : s , s, 1) ; 

A=l/s*(j<i)  +  l/(2*s)*(i==j) ; 
b=l/s*ones(s, 1) ; 
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c=sum(A,2) ; 


case  ’ISSPs3’ 

‘/.Optimal  DIRK  SSP  schemes  of  order  3 

if  s<2  error (’ Implicit  third  order  SSP  schemes  require  s>=2’);  end 
r=s-l+sqrt(s"2-l) ; 

i=repmat ( (1 : s) ’ , 1 , s) ;  j=repmat (1 ; s , s , 1) ; 

A=l/sqrt(s"2-l)*(j<i)  +  0 . 5* (1-sqrt ( (s-1) /(s+1) ) ) * (i==j ) ; 
b=l/s*ones(s, 1) ; 
c=sum(A,2) ; 


Classical  Methods' 


‘/.Gauss-Legendre  methods  —  order  2s 
case  'GLl' 

r=2;  A=l/2;  b=l;  c=l/2; 
case  ’GL2’ 
r=0; 

A=[l/4  l/4-sqrt(3)/6 
l/4+sqrt(3)/6  1/4]; 
b=[l/2  1/2]’; 

c= [l/2-sqrt(3)/6  l/2+sqrt(3)/6] ’ ; 
case  ’GL3’ 
r=0; 

A=[5/36  (80-24*sqrt(15))/360  (50-12*sqrt(15) )/360 
(50+15*sqrt(15))/360  2/9  (50-15*sqrt(15) )/360 

(50+12*sqrt(15))/360  (80+24*sqrt (15) ) /360  5/36]; 
b=[5/18  4/9  5/18] ’ ; 

c=[(5-sqrt(15))/10  1/2  (5+sqrt (15) ) /lO] ’ ; 

'/.Radau  lA  methods  —  order  2s-l 
case  ’RIAl’ 
r=l; 

A=l;  b=l;  c=0; 
case  ’RIA2’ 
r=0; 

A=[l/4  -1/4 
1/4  5/12] ; 
b=[l/4  3/4] ’ ; 
c=[0  2/3]  ’  ; 
case  ’RIA3’ 
r=0; 

A=[l/9  (-l-sqrt(6))/18  (-l+sqrt(6))/18 

1/9  (88+7*sqrt(6))/360  (88-43*sqrt(6))/360 
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1/9  (88+43*sqrt(6))/360  (88-7*sqrt(6) )/360] ; 
b=[l/9  (16+sqrt(6))/36  (16-sqrt(6) )/36] ’ ; 
c=[0  (6-sqrt(6) )/10  (6+sqrt(6))/10] ' ; 

’/.Radau  IIA  methods  —  order  2s- 1 
case  ’RIIAl’ 
r=l; 

A=l;  b=l;  c=l; 
case  'RIIA2’ 
r=0; 

A=[5/12  -1/12 
3/4  1/4]; 
b=[3/4  1/4]’; 
c=[l/3  1]  ’  ; 
case  ’RIIA3’ 
r=0; 

A=[(88-7*sqrt(6))/360  (296-169*sqrt(6))/1800  (-2+3*sqrt(6))/225 
(296+169*sqrt(6))/1800  (88+7*sqrt (6) ) /360  (-2-3*sqrt(6))/225 
(16-sqrt(6))/36  (16+sqrt(6) )/36  1/9]; 
b=[(16-sqrt(6))/36  (16+sqrt(6))/36  1/9  ]’; 
c=[(4-sqrt(6))/10  (4+sqrt(6))/10  1]; 


y.Lobatto  I  IIA  methods  —  order  2s-2 
case  ’LIIIA2’ 
r=0; 

A=[0  0 

1/2  1/2] ; 
b=[l/2  1/2] ’ ; 
c=[0  1]’; 
case  ’LIIIA3’ 
r=0; 

A=[0  0  0 

5/24  1/3  -1/24 
1/6  2/3  1/6] ; 
b=[l/6  2/3  1/6] ’ ; 
c= [0  12  1]  ; 


Miscellaneous  Methods 


case  ’Mid22’ 

"/.Midpoint  22  method 
s=2;  r=0.5; 

A=[0  0 

1/2  0] ; 

b=[0  1]  ’  ;  c=[0  1/2]  ’  ; 
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case  'MTE22' 

*/, Minimal  truncation  error  22  method  (Heun) 
s=2 ;  r=0 . 5 ; 

A=[0  0 

2/3  0]; 

b=[l/4  3/4] ' ;  c=[0  2/3] ' ; 

case  'CN22' 

y.Crank-Nicholson 
s=2 ;  r=2 ; 

A=[0  0 

1/2  1/2]; 

b=[l/2  1/2]  '  ;  c=[0  1]  '  ; 

case  ’Heun33’ 
s=3;  r=0; 

A=[0  0  0;  1/3  0  0;  0  2/3  0]; 
b=[l/4  0  3/4]’;  c=sum(A,2); 

case  ’RK44’  '/.Classical  fourth  order 
s=4 ;  r=0 ; 

A=[0  0  0  0;  1/2  0  0  0;  0  1/2  0  0;  0  0  1  0] 
b=[l/6  1/3  1/3  1/6]’;  c=sum(A,2); 

*/.====================DSRK  Methods============= 


case  ’DSso2’ 

•/.CBM’s  DSRKso2 
s=2;  isdsrk=l; 

A=[3/4  -1/4 

1  0]  ; 

W=[l/2  0 

1  0]  ; 

b=[l  0]  ’  ;  c=[l/2  1]  ’  ; 

case  ’DSRK2’ 

•/.CBM’s  DSRK2 
s=2;  isdsrk=l; 

A=[l/2  -1/2 

1/2  1/2] ; 

W=[  0  0 

1/2  1/2] ; 

b=[l/2  1/2]’;  c=[0  1]’; 
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case  'DSRK3' 

'/.Zennaro's  DSRK3 
s=3;  isdsrk=l; 

A=[5/2  -2  -1/2 
-1  2  -1/2 
1/6  2/3  1/6] ; 

W=[  0  0  0 

7/24  1/6  1/24 
1/6  2/3  1/6] ; 

b=[l/6  2/3  1/6]';  c=[0  1/2  1]'; 

’/,==================== "Non-SSP"  Methods  of  Wong  &  Spiteri 

case  'NSSP21' 
m=2 ;  r=0 ; 

A=[0  0 

3/4  0] ; 

b=[0  1]  '  ;  c=[0  3/4]  '  ; 

case  'NSSP32’ 
in=3;  r=0 ; 

A=  [0  0  0 

1/3  0  0 
0  10]; 

b=[l/2  0  1/2]  '  ;  c=[0  1/3  1]  '  ; 

case  'NSSP33' 
m=3;  r=0; 

A=[0  0  0 

-4/9  0  0 

7/6  -1/2  0] ; 

b=[l/4  0  3/4]';  c=[0  -4/9  2/3]'; 

case  'NSSP53' 
m=5;  r=0; 

A=[0  0  0  0  0 
1/70000 
0  3/16  000 
001/300 
0  0  0  2/3  0]; 

b=[l/4  000  3/4]';  c=[0  1/7  3/16  1/3  2/3]'; 

end 


•  Order  of  a  Runge  -Kutta  method: 


function  p=rk_order(A,b,c) 
’/.By  David  Ketcheson 
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7o 

7oDetermine  order  of  a  RK  method,  up  to  sixth  order 
7oOrder  conditions  from  text  of  Hairer,  Norsett,  &  Wanner 

7« 

"/oFor  an  m-stage  method,  input  A  should  be  a  m  x  m  matrix; 
"/ob  and  c  should  be  column  vectors  of  length  m 

eps=l . e-14; 

m=length(b);  7o  #  of  stages 

em=ones(m, 1) ; 

p=0; 

if  sum(b)“l<eps 

p=i; 

end 

z(l)=sum(A^  *b)-l/2; 

if  (p==l  &&  abs(z(l))<eps)  p=2;  end 

z(l)=c^ .^2*b-l/3; 

z(2)=b'  *A''2*em”l/6 ; 

if (max(abs(z))<eps  &&  p==2)  p=3;  end 

z(l)=b'*c. ^3-1/4; 

z(2)=(b . *c) ' *A"2*ones(m, l)-l/8; 

z(3)=b^*A*c. '‘2-1/12; 

z(4)=b^*A'‘2*c-l/24; 

if (max(abs(z))<eps  &&  p==3)  p=4;  end 

z(l)=c^^4*b-l/5; 
z(2)=(b.*c.^2) '*A*c-l/10; 
z(3)=b'*(A*c)  .'‘2-1/20; 
z(4)=(b . *c) ’ *A*c . "2-1/15 ; 

2(5)=b'*A*c. "3-1/20; 

z(6)=(b.*c)^*A"2*c-l/30; 

z(7)=b ’ *A*diag(c) *A*c-l/40; 

z(8)=b^*A"2*c. "2-1/60; 

z(9)=b^*A"3*c-l/120; 

if (max(abs(z))<eps  &&  p==4)  p=5;  end 

if  p==5 

z(l)=c^"5*b-l/6; 
z(2)=b^*diag(c) . "3*A*c-l/12; 
z(3)=b^*diag(c)*(A*c) . "2-1/24; 
z(4)=b^ *diag(c) . "2*A*c. "2-1/18; 
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z(5)=b^*((A*c.^2) .*(A*c))-l/36; 
z(6)=b^  *diag(c)*A*c . "3-1/24; 
z(7)=b^*A*c. "4-1/30; 
z(8)=b^*diag(c) . "2*A"2*c-l/36; 
z(9)=b^*((A"2*c) .*(A*c))-l/72; 
z(10)=b^*diag(c)*A*diag(c)*A*c-l/48; 
z(ll)=b^*A*diag(c) . "2*A*c-l/60; 
z(12)=b^*A*(A*c) ."2-1/120; 
z(13)=b^  *diag(c)*A"2*c . "2-1/72; 
z(14)=b^  *A*diag(c) *A*c . "2-1/90; 
z(15)=b^*A"2*c. "3-1/120; 
z(16)=b' *diag(c) *A"3*c- 1/144; 
z(17)=b'*A*diag(c)*A"2*c-l/180; 
z(18)=b^  *A"2*diag(c) *A*c-l/240 ; 
z(19)=b^*A"3*c. "2-1/360; 
z(20)=b^*A"4*c-l/720; 

if (max(abs(z))<eps)  p=6;  print(^This  method  has  order  at  least  sixO;  end 
end 

•  Absolute  monotonic  polynomials: 

f unct ion  [gamma , R] =Rsp ( s , p) 

7oBy  David  Ketcheson 

y. 

"/oReturns  the  optimal  absolutely  monotonic  polynomial  of  degree  s 
"/and  order  of  accuracy  p 

"/gamma  contains  the  coefficients  of  the  Taylor  series  about  z=-r 
"/To  construct  the  polynomial,  use: 

"/>  syms  z  phi 

"/>  phi=simplify  (sum(  (1 .  +z/R)  .  "  (0  :  s)  .  *gamma) ) ; 

“/ 

"/Uses  the  MATLAB  optimization  toolbox 


%========================================================= 

"/Set  options  for  linprog 

opts=optimset ( 'TolX^ , 1 . e-15 , ^TolFim^ , 1 . e-15 , ^Maxlter ' , 10000000 , . . . 

^LargeScale' , ^on\  ^ Simplex^ , ^off ^ ^ Display^ , ^off 0 ; 
acc=l.e-15;  "/Accuracy  of  bisection  search 
%========================================================= 


if  p==s  "/In  this  case,  the  optimal  polynomial  is  just  the  Taylor  polynomial 
R=l; 

for  i=0:p 
d(i+l)=R"i ; 
for  j=0:s 
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B(i+1, j+l)=prod(j-(0: i-1)) ; 
end 
end 

gainma=(B\ones(s+l ,  1))  ^  ; 

else 
M=s+1 ; 

rmax=S“'p+l .  0001 ; 
rmin=0; 

r=rmax;  "/plnitial  guess 
c=zeros(M, 1) ; 
clear  B  d; 

while  (rmax-rmin>acc)  "/oFind  R  by  bisection 

'/oSet  up  and  improve  conditioning  of  equality  constraints 
for  i=0:p 

rescale=r''i ;  d(i+l)=r"i/rescale ; 
for  j=0:s 

B(i+1 , j+l)=prod( j-(0 : i-1) ) /rescale ; 
end 
end 

"/oTest  feasibility  for  this  value  of  r 

[x, lambda, exitf lag] =linprog(c , []  , []  ,B,d,zeros(M, 1) ,zeros(M, 1)+1 .e6,c,opts) ; 
if  exitf lag==l; 

rmin=r ;  r= (r+rmax) /2 ; 
else 

rmax=r ;  r= (rmin+r ) / 2 ; 
end 
end 

"/pNow  get  a  feasible  solution  so  we  have  the  coefficients  of  the  method 
R=rmin; 
for  i=0:p 
rescale=R''i ; 
d(i+l)=R''i/rescale ; 
for  j=0:s 

B(i+1 , j+l)=prod(j-(0 : i-1) ) /rescale ; 
end 
end 

[gamma, lambda, exit f lag] =linprog(c, []  , []  , B,d, zeros (M, 1) , zeros (M, 1)  +  1 . e6 , c , opts) 
end 
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